#include "annealing.h"

#ifdef ENABLE_MULTIMAX

t_annealing anneal;

// "Set_Anneal" sets the default values of thermal annealing parameters. . .
void Set_Annealing_Defaults(option *io){
	anneal.t_scalar = 0.3; //this one gets set from io later... don't change this line.
	anneal.random_seed = (int)time(NULL);
	PhyML_Printf("\n. Simulated annealing random seed = %d\n", anneal.random_seed);
	anneal.rng = gsl_rng_alloc(gsl_rng_mt19937);
	anneal.no_change = 1;
	anneal.max_alpha = 0.3;
	gsl_rng_set(anneal.rng, anneal.random_seed);
}

void Free_Anneal(){
	gsl_rng_free( anneal.rng );
}


void Swap_Tree_And_Mod(t_tree *tree1, t_tree *tree2){
	t_tree *tmptree;
	model *tmpmod;
	tmptree = tree1;
	tmpmod = tree1->mod;

	tree1 = tree2;
	tree1->mod = tree2->mod;
	tree2 = tmptree;
	tree2->mod = tree1->mod;
}

// This method returns the scaled and estimated acceptance ratio for thermal annealing.
// The variable "tree" will be modified by this method, such that tree will be in a good
// starting position for thermal annealing after this method concludes.
phydbl Scale_T_Factor(t_tree *tree){
	if (!tree->io->quiet)
	{
		PhyML_Printf("\n. Scaling temperature.  This could take a while...\n");
	}
	phydbl aratio;
	Lk(tree);
	phydbl lnL_best = tree->c_lnL;
	phydbl lnL_proposed = tree->c_lnL;
	phydbl fsum = 0.0;
	phydbl fsqsum = 0.0;
	int n = 0;
	int i;

	t_tree *best_tree = best_tree = Make_Tree_From_Scratch(tree->n_otu,tree->data);
	best_tree->mod  = Copy_Model(tree->mod);
#ifndef ENABLE_MBL
	phydbl *best_bl = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
#else
	phydbl *best_bl = (phydbl *)mCalloc((2*tree->n_otu-3)*n_l,sizeof(phydbl));
#endif
	Record_Br_Len(best_bl,tree);
	Copy_Tree(tree, best_tree);

	For(i, tree->io->sa_iters_per_stage){
		Get_Neighbor_Proposition(tree,tree->io->sa_temp_start);
		lnL_proposed = tree->c_lnL;
		if(lnL_proposed > UNLIKELY){ // is the proposition valid?
			n++;
			fsum += lnL_proposed;
			fsqsum += (lnL_proposed * lnL_proposed);
		}
		if(lnL_proposed > lnL_best){ // is the proposition better?
			lnL_best = lnL_proposed;
			Copy_Tree(tree,best_tree); // write tree into best_tree
			Record_Model(tree->mod, best_tree->mod); // write tree->mod into best_tree->mod
			Record_Br_Len(best_bl, tree);
		}else{
			Copy_Tree(best_tree,tree);
			Record_Model(best_tree->mod, tree->mod);
			Restore_Br_Len(best_bl, tree);
		}

		if (!tree->io->quiet)
		{
			phydbl a = fsum / (phydbl) n;
			phydbl b = fsqsum / (phydbl) n;
			phydbl c = sqrt(b - a * a);
			PhyML_Printf(" . [ iteration %d or %d] [scalar = %f]\n", i, tree->io->sa_iters_per_stage, anneal.t_scalar * c / tree->io->sa_temp_start);
		}
	}
	fsum /= (phydbl)n; //average proposed likelihood
	fsqsum /= (phydbl)n; //average (squared) proposed likelihood
	fsqsum = sqrt(fsqsum - fsum * fsum);
	aratio = anneal.t_scalar * fsqsum / tree->io->sa_temp_start;

	Free_Model(best_tree->mod);
	Free(best_bl);
	Free_Tree(best_tree);

	//now lets get a really good starting position...
	//Speed_Spr_Loop(tree);

	return aratio;
}

// This method modifies some or all parameters of "tree", including:
// * topology
// * branch lengths
// * branch length proportions
// * alpha (for gamma distributed ASRV)
// * gamma proportions
// * mu (evolutionary rate)
void Get_Neighbor_Proposition(t_tree *tree,phydbl temp){

#ifdef TIME
	clock_t time_STA_GetNeighborProposition_start = clock();
#endif

	tree->both_sides = 0;
	anneal.no_change = 1;
	phydbl x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_topo)Step_Topology(tree);
	x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_kappa)Step_Kappa(tree);
	x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_lambda)Step_Lambda(tree);
	x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_rr)Step_RR(tree);
	x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_pi) Step_Pi(tree);
	x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_gamma)Step_Gamma(tree);
	x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_pinvar)Step_Pinvar(tree);
	x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_brlen)Step_Branch_Lengths(tree,temp);
#ifdef ENABLE_MBL
	x = gsl_rng_uniform(anneal.rng);
	if(x < tree->io->sa_prob_blprop)Step_Brlen_Proportion(tree);
#endif

	// Update the likelihood of tree, but only if the model was changed.
	if(anneal.no_change == 0)
	{
		Lk(tree);
	}
	tree->mod->update_eigen = 0;
	tree->both_sides = 1; // reset to 1

	//
	// Print the full proposition to a text file
	//
	if (tree->io->print_trace == 1)
	{
#ifndef ENABLE_MBL
		char *s = Write_Tree(tree);
#else
		char *s = Write_Tree(tree,-1);
#endif
		PhyML_Fprintf(tree->io->fp_out_trace,"[%f]%s\n",tree->c_lnL,s); fflush(tree->io->fp_out_trace);
		Free(s);

		int i;
		PhyML_Fprintf(tree->io->fp_out_annealing_file,
				". proposition: [lnL = %f]\n",
				tree->c_lnL);
#ifdef ENABLE_MBL
		// BL proportions
		if(tree->mod->s_opt->opt_bl_mixtures){
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "\t[BL proportions] [");
			For(i, n_l)
			{
				PhyML_Fprintf(tree->io->fp_out_annealing_file, "%.3f ", tree->mod->bl_mixture_wts[i]);
			}
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "]\n");
		}
#endif

		// tree->mod->alpha
		if(tree->mod->s_opt->opt_alpha == 1){
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "\t[alpha] [%.3f]\n", tree->mod->alpha);
		}

		// tree->mod->pinvar
		if(tree->mod->s_opt->opt_pinvar == 1){
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "\t[prob. invar.] [%.3f]\n", tree->mod->pinvar);
		}

		// tree->mod->kappa
		if(tree->mod->s_opt->opt_kappa == 1){
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "\t[kappa] [%.3f]\n", tree->mod->kappa);
		}

		// tree->mod->lambda
		if(tree->mod->s_opt->opt_lambda == 1){
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "\t[lambda] [%.3f]\n", tree->mod->lambda);
		}

		// relative rate params. for GTR
		if(((tree->mod->whichmodel == GTR) && (tree->mod->s_opt->opt_rr == 1)) ||
				((tree->mod->whichmodel == CUSTOM) && (tree->mod->s_opt->opt_rr == 1)) )
		{
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "\t[relative rates] [");
			for(i = 0; i < tree->mod->n_diff_rr; i++){
				PhyML_Fprintf(tree->io->fp_out_annealing_file, "%.3f ", tree->mod->rr_val[i]);
			}
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "]\n");
		}

		// state frequencies
		if(tree->mod->s_opt->opt_state_freq)
		{
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "\t[state freqs.] [");
			for(i = 0; i < tree->mod->ns; i++){
				PhyML_Fprintf(tree->io->fp_out_annealing_file, "%.3f ", tree->mod->pi_unscaled[i]);
			}
			PhyML_Fprintf(tree->io->fp_out_annealing_file, "]\n");
		}

		fflush(tree->io->fp_out_annealing_file);
	}

#ifdef TIME
	clock_t time_STA_GetNeighborProposition_stop = clock();
	time_STA_GetNeighborProposition += time_STA_GetNeighborProposition_stop - time_STA_GetNeighborProposition_start;
#endif
}


#ifdef ENABLE_MBL
void Step_Brlen_Proportion(t_tree *tree){
	if(tree->mod->s_opt->opt_bl_mixtures)
	{
		if (tree->io->sa_sele_blprop == 1)
		{
			int max = tree->io->sa_sele_max_params;
			if (max < 1) max = 1;
			int n_to_modify = gsl_rng_uniform_int(anneal.rng,max);
			if(n_to_modify < 1) n_to_modify = 1; // we'll perturb the values of this many BL props

			int i;
			For(i,n_to_modify)
			{
				// select a random BL proportion value (r_i)
				int r_i = gsl_rng_uniform_int(anneal.rng, n_l);

				phydbl r;
				// select a random perturbation value, r
				//r = gsl_ran_gaussian(anneal.rng,tree->io->sa_sigma_blprop);
				r = gsl_rng_uniform(anneal.rng) * tree->io->sa_sigma_blprop;
				int p = gsl_rng_uniform_int(anneal.rng, 2); // to make rand_gauss + or -?
				if ( p == 1 ) r = r * -1;
				tree->mod->bl_mixture_wts[r_i] += r;
				if(tree->mod->bl_mixture_wts[r_i] < 0.001) tree->mod->bl_mixture_wts[r_i] = 0.001;
				else if(tree->mod->bl_mixture_wts[r_i] > 100.0) tree->mod->bl_mixture_wts[r_i] = 99.9999;
			}
		}
		else if (tree->io->sa_sele_blprop == 2)
		{
            phydbl *alpha;
            alpha = (phydbl *)mCalloc(n_l, sizeof(phydbl));
            int i;
            For(i,n_l){
            	alpha[i] = anneal.max_alpha;
            }
            gsl_ran_dirichlet(anneal.rng,n_l,alpha,tree->mod->bl_mixture_wts);
            Free(alpha);
		}

		// and finally...
        tree->both_sides = 1;
        anneal.no_change = 0;
        Set_nth_mixture_weight(tree->mod->bl_mixture_wts, n_l);
	}
}
#endif


void Step_Pinvar(t_tree *tree){
	if(tree->mod->s_opt->opt_pinvar == 1){
		if (tree->io->sa_sele_pinvar == 0) // Step using a hill-climbing
		{
			Optimize_Single_Param_Generic(tree,&(tree->mod->pinvar),
					.0001,0.9999,
					tree->mod->s_opt->min_diff_lk_global,
					tree->mod->s_opt->brent_it_max,
					tree->mod->s_opt->quickdirty);
		}
		if (tree->io->sa_sele_pinvar == 1)
		{
			phydbl r;
			anneal.no_change = 0;
			//r = gsl_ran_gaussian(anneal.rng,tree->io->sa_sigma_pinvar);
			r = gsl_rng_uniform(anneal.rng) * tree->io->sa_sigma_pinvar;
			int p = gsl_rng_uniform_int(anneal.rng, 2); // to make rand_gauss + or -?
			if ( p == 1 ) r = r * -1;
			tree->mod->pinvar += r;
			if(tree->mod->pinvar < 0.0001) tree->mod->pinvar = 0.0001;
			else if(tree->mod->pinvar > 0.9999) tree->mod->pinvar = 0.9999;
		}
	}
}

void Step_RR(t_tree * tree){
	if(((tree->mod->whichmodel == GTR) && (tree->mod->s_opt->opt_rr == 1) ) ||
			((tree->mod->whichmodel == CUSTOM) && (tree->mod->s_opt->opt_rr == 1)))
	{
		if (tree->io->sa_sele_rr == 0) // Step using a hill-climbing
		{
			  int i;
			  For(i,tree->mod->n_diff_rr)
			  {
				  if(i != 5)
				  {
					  Generic_Brent_Lk(&(tree->mod->rr_val[i]),
						 1.E-2,1.E+2,
						 tree->mod->s_opt->min_diff_lk_global,
						 tree->mod->s_opt->brent_it_max,
						 tree->mod->s_opt->quickdirty,
						 Wrap_Lk,NULL,tree,NULL);
				  }
			  }
		}
		else if (tree->io->sa_sele_rr == 1)
		{
			int r_i = gsl_rng_uniform_int(anneal.rng, tree->mod->n_diff_rr);// select a random RR value to perturb
			phydbl r;
			anneal.no_change = 0;
			tree->mod->update_eigen = 1;
			//r = gsl_ran_gaussian(anneal.rng,tree->io->sa_sigma_rr);
			r = gsl_rng_uniform(anneal.rng) * tree->io->sa_sigma_rr;
			int p = gsl_rng_uniform_int(anneal.rng, 2); // to make rand_gauss + or -?
			if ( p == 1 ) r = r * -1;
			tree->mod->rr_val[r_i] += r;
			if(tree->mod->rr_val[r_i] < 0.001) tree->mod->rr_val[r_i] = 0.001;
			else if(tree->mod->rr_val[r_i] > 100.0) tree->mod->rr_val[r_i] = 99.9999;
		}

	}
}


void Step_Kappa(t_tree * tree){
	if(tree->mod->s_opt->opt_kappa == 1){
		if (tree->io->sa_sele_kappa == 0) // Step using a hill-climbing
		{
			tree->mod->update_eigen = 1;
			anneal.no_change = 0;
			Optimize_Single_Param_Generic(tree,&(tree->mod->kappa),
					.1,100.,
					tree->mod->s_opt->min_diff_lk_global,
					tree->mod->s_opt->brent_it_max,
					tree->mod->s_opt->quickdirty);
		}
		else if (tree->io->sa_sele_kappa == 1)
		{
			phydbl r;
			tree->mod->update_eigen = 1;
			anneal.no_change = 0;
			//r = gsl_ran_gaussian(anneal.rng,tree->io->sa_sigma_kappa);
			r = gsl_rng_uniform(anneal.rng) * tree->io->sa_sigma_kappa;
			int p = gsl_rng_uniform_int(anneal.rng, 2); // to make rand_gauss + or -?
			if ( p == 1 ) r = r * -1;
			tree->mod->kappa += r;
			if(tree->mod->kappa < 0.1) tree->mod->kappa = 0.1;
			else if(tree->mod->kappa > 100.0) tree->mod->kappa = 99.9999;
		}
	}
}


void Step_Lambda(t_tree * tree){
	if(tree->mod->s_opt->opt_lambda == 1){
		if (tree->io->sa_sele_lambda == 0) // Step using a hill-climbing
		{
			anneal.no_change = 0;
			tree->mod->update_eigen = 1;
			Optimize_Single_Param_Generic(tree,&(tree->mod->lambda),.001,100.,
					tree->mod->s_opt->min_diff_lk_global,
					tree->mod->s_opt->brent_it_max,
					tree->mod->s_opt->quickdirty);
		}
		else if (tree->io->sa_sele_lambda == 1)
		{
			phydbl r;
			anneal.no_change = 0;
			tree->mod->update_eigen = 1;
			//r = gsl_ran_gaussian(anneal.rng,tree->io->sa_sigma_lambda);
			r = gsl_rng_uniform(anneal.rng) * tree->io->sa_sigma_lambda;
			int p = gsl_rng_uniform_int(anneal.rng, 2); // to make rand_gauss + or -?
			if ( p == 1 ) r = r * -1;
			tree->mod->lambda += r;
			if(tree->mod->lambda < 0.001) tree->mod->lambda = 0.001;
			else if(tree->mod->lambda > 100.0) tree->mod->lambda = 99.9999;
		}
	}
}


void Step_Gamma(t_tree *tree){
	if((tree->mod->s_opt->opt_alpha == 1) && (tree->mod->n_catg > 1)){
		if (tree->io->sa_sele_gamma == 0) // Step using a hill-climbing
		{
			anneal.no_change = 0;
			Optimize_Single_Param_Generic(tree,&(tree->mod->alpha),
					tree->mod->alpha/2.,tree->mod->alpha*2.,
					tree->mod->s_opt->min_diff_lk_global,
					tree->mod->s_opt->brent_it_max,
					tree->mod->s_opt->quickdirty);
		}
		else if (tree->io->sa_sele_gamma == 1)
		{
			anneal.no_change = 0;
			//phydbl r = gsl_ran_gaussian(anneal.rng,tree->io->sa_sigma_gamma);
			phydbl r = gsl_rng_uniform(anneal.rng) * tree->io->sa_sigma_gamma;
			int p = gsl_rng_uniform_int(anneal.rng, 2); // to make rand_gauss + or -?
			if ( p == 1 ) r = r * -1;
			tree->mod->alpha += r;
			if(tree->mod->alpha < 0.01) tree->mod->alpha = 0.01;
			else if(tree->mod->alpha > 100.0) tree->mod->alpha = 100.0;
		}
	}
}

void Step_Pi(t_tree * tree){
	if((tree->mod->s_opt->opt_state_freq))
	{
		if (tree->io->sa_sele_pi == 1)
		{
			int max = tree->io->sa_sele_max_params;
			if (max < 1) max = 1;
			int n_pi_to_modify = gsl_rng_uniform_int(anneal.rng,max);
			if(n_pi_to_modify < 1) n_pi_to_modify = 1; // we'll perturb the values of this many pi vals

			int i;
			For(i,n_pi_to_modify){
				// select a random pi value (r_i)
				phydbl x = gsl_rng_uniform(anneal.rng);
				int r_i = (int)x%tree->mod->ns;

				phydbl r;
				anneal.no_change = 0;
				tree->mod->update_eigen = 1;
				// select a random perturbation value, r
				//r = gsl_ran_gaussian(anneal.rng,tree->io->sa_sigma_pi);
				r = gsl_rng_uniform(anneal.rng) * tree->io->sa_sigma_pi;
				int p = gsl_rng_uniform_int(anneal.rng, 2); // to make rand_gauss + or -?
				if ( p == 1 ) r = r * -1;
				//PhyML_Printf(". debug: annealing.c:452, proposing new gaussian, r = %f.\n", r);
				tree->mod->pi_unscaled[r_i] += r;
				if(tree->mod->pi_unscaled[r_i] < 0.001) tree->mod->pi_unscaled[r_i] = 0.001;
				else if(tree->mod->pi_unscaled[r_i] > 100.0) tree->mod->pi_unscaled[r_i] = 99.9999;
			}
		}
		if (tree->io->sa_sele_pi == 0) // Step using a hill-climbing
		{
			int i;
			For(i, tree->mod->ns)
		    {
		    	Generic_Brent_Lk(&(tree->mod->pi_unscaled[i]),
					 -1000.,1000.,
					 tree->mod->s_opt->min_diff_lk_global,
					 tree->mod->s_opt->brent_it_max,
					 tree->mod->s_opt->quickdirty,
					 Wrap_Lk,NULL,tree,NULL);
		    }
			tree->mod->update_eigen = 0;
		}
		else if (tree->io->sa_sele_pi == 2) // Step using a Dirichlet distribution
		{
			tree->mod->update_eigen = 1;
			anneal.no_change = 0;
            phydbl *alpha;
            alpha = (phydbl *)mCalloc(tree->mod->ns, sizeof(phydbl));
            int i;
            For(i,tree->mod->ns){
            	alpha[i] = anneal.max_alpha;
            }
            gsl_ran_dirichlet(anneal.rng,tree->mod->ns,alpha,tree->mod->pi_unscaled);
            Free(alpha);
		}
	}
}

// helper for "Get_Neighbor_Proposition"
void Step_Branch_Lengths(t_tree *tree, phydbl temp){
	if(tree->mod->s_opt->opt_bl == 1)
	{
		if (tree->io->sa_sele_bl == 0) // Step using a hill-climbing
		{
			Optimize_Br_Len_Serie(tree->noeud[0],
									   tree->noeud[0]->v[0],
									   tree->noeud[0]->b[0],
									   tree,
									   tree->data);
			tree->both_sides = 1;
			anneal.no_change = 0;
		}
		else if(tree->io->sa_sele_bl == 1)
		{
			tree->both_sides = 1;
			anneal.no_change = 0;
			int i,j,n_edges_to_modify,n,max;
			int n_edges = (tree->n_otu * 2) - 3;
			phydbl rand_gauss;

			max = tree->io->sa_sele_max_params;
			if (max < 1) max = 1;
			n_edges_to_modify = gsl_rng_uniform_int(anneal.rng,max);
			if(n_edges_to_modify < 1) n_edges_to_modify = 1; // we'll perturb the lengths of this many edges

			For(i,n_edges_to_modify){
				j = gsl_rng_uniform_int(anneal.rng,n_edges); // j is a random edge


				// Victor's note: I'm having difficulty getting gsl_ran_gassuain to work correctly.
				// see: http://flamingovic.wordpress.com/2011/03/17/random-gaussian-distribution-in-libgsl/
				//rand_gauss = gsl_ran_gaussian(anneal.rng, (double)tree->io->sa_sigma_bl); // get a random perturbation delta
				// This is a hack for now:
				rand_gauss = gsl_rng_uniform(anneal.rng) * tree->io->sa_sigma_bl;
				int p = gsl_rng_uniform_int(anneal.rng, 2); // to make rand_gauss + or -?
				if ( p == 1 ) rand_gauss = rand_gauss * -1;
#ifdef ENABLE_MBL
				n = gsl_rng_uniform_int(anneal.rng,n_l); // n is a random BL set
				tree->t_edges[j]->l[n] += rand_gauss; // perturb branch j's nth value
				//PhyML_Printf(" . debug annealing.c edge %d, mixture %d, length += %f\n", j, n, rand_gauss);
				if(tree->t_edges[j]->l[n] < BL_MIN) tree->t_edges[j]->l[n] = BL_MIN;
				else if(tree->t_edges[j]->l[n] > BL_MAX) tree->t_edges[j]->l[n] = BL_MAX;
#else
				tree->t_edges[j]->l += rand_gauss; // perturb branch j's nth value
				//PhyML_Printf(" . debug annealing.c 513 edge %d, sigma = %f, BL stepsize = %f\n", j, tree->io->sa_sigma_bl, rand_gauss);
				if(tree->t_edges[j]->l < BL_MIN) tree->t_edges[j]->l = BL_MIN;
				else if(tree->t_edges[j]->l > BL_MAX) tree->t_edges[j]->l = BL_MAX;
#endif

			}
		}
	}
}


// helper for "Get_Neighbor_Proposition"
void Step_Topology(t_tree *tree){
	if(tree->mod->s_opt->opt_topo)
	{
		// Make nmoves number of topology steps, using either NNI or SPR
		int nmoves = gsl_rng_uniform_int(anneal.rng, 10);

		//phydbl p = gsl_rng_uniform(anneal.rng);
		//if(p <= tree->io->sa_prob_NNI)
		//{
		  Random_NNI(nmoves,tree);
		  anneal.no_change = 0;
		//}

		//p = gsl_rng_uniform(anneal.rng);
		//if(p <= tree->io->sa_prob_SPR)
		//{
		//	Random_Spr(nmoves,tree);
		//	anneal.no_change = 0;
		//}
	}
}


void Get_QA_Neighbor_Proposition(t_tree *tree){

}

// Boltzmann probability sampler, for thermal annealing
//
// INPUT: Given the likelihood (lnl_curr) of our current state, the likelihood (lnl_new) of our new proposed state,
// and the temperature (temperature) of the system, calculate the probability of the new state,
// where the probability is drawn from a Boltzmann distribution.
//
// OUTPUT: a floating point decimal, between 0.0 and 1.0
phydbl Boltzmann_P_TA(phydbl lnl_curr, phydbl lnl_new, phydbl temperature){
	phydbl result = 1.0;
	if(lnl_new < lnl_curr){
		result = exp((lnl_new - lnl_curr)/(anneal.t_scalar * temperature));
	}
	//PhyML_Printf("In Boltzmann_P: t_scalar= %lf, lnl_curr = %lf, lnl_new = %lf, temp= %lf, result = %lf\n",anneal.t_scalar,lnl_curr,lnl_new,temperature,result);
	return result;
}


// INPUT: a tree structure, where tree->mod->s_opt contains the
//      parameters which are free to be optimized.
// OUTPUT: the likelihood of the best found tree
//
// At the end of this method, the object "tree" will contain the best-found topology,
//      branch lengths, and model parameters.  In other words, tree is mutable!
phydbl Thermal_Anneal_All_Free_Params(t_tree *tree, int verbose){

#ifdef TIME
	clock_t time_STAstart = clock();
#endif

	Set_Annealing_Defaults(tree->io);
	tree->mod->update_eigen = 1;
	tree->both_sides = 1;
	int n_edges = (tree->n_otu * 2) - 3;
	if (n_edges <= 3){
		tree->mod->s_opt->opt_topo = 0; //make sure that opt_topo is false if there are no meaningful branch swaps.
	}

	//
	// Pre-optimize the tree?
	//
	if (tree->io->sa_preopt > -1)
	{
		if (!tree->io->quiet)
		{	PhyML_Printf("\n. Before STA begins: pre-optimizing topology, branch lengths, and model parameters.\n");
		}

		tree->io->opt_algorithm = tree->io->sa_preopt;
		if(tree->mod->s_opt->opt_topo)
		{
			  if(tree->mod->s_opt->topo_search      == NNI_MOVE) Simu_Loop(tree);
			  else if(tree->mod->s_opt->topo_search == SPR_MOVE) Speed_Spr_Loop(tree);
			  else                                               Best_Of_NNI_And_SPR(tree);
		}
		else
		{
			  if(tree->mod->s_opt->opt_subst_param ||
			 tree->mod->s_opt->opt_bl)                       Round_Optimize(tree,tree->data,ROUND_MAX);
			  else                                               Lk(tree);
		}

		//
		// Print the pre-opt tree:
		//
		char pfilestr[1000];
		sprintf(pfilestr, "%s.preopt", tree->io->out_tree_file);
		FILE* pfile = fopen(pfilestr, "w");
		if(!pfile)
		{
			fprintf(stderr, "ERROR: cannot open output file: %s\n",
						pfilestr);
			exit(1);
		}
#ifndef ENABLE_MBL
		char *most_likely_trees = Write_Tree(tree);
#else
		char *most_likely_trees = Write_Tree(tree, -1);
#endif
		PhyML_Fprintf(pfile,"%s\n",most_likely_trees);
		fclose(pfile);
		Free(most_likely_trees);

		//
		// Print the pre-opt stats file:
		//
		char pfilestr2[1000];
		sprintf(pfilestr2, "%s.preopt", tree->io->out_stats_file);
		pfile = fopen(pfilestr2, "w");
		if(!pfile)
		{
			fprintf(stderr, "ERROR: cannot open output file: %s\n",
						pfilestr2);
			exit(1);
		}
		time_t t_end;
		time(&t_end);
		Print_Fp_Out(pfile,tree->t_beg,t_end,tree,tree->io,1,1);
		fclose(pfile);

		//
		// Print the pre-optimized per-site lnls
		//
		char pfilestr3[1000];
		sprintf(pfilestr3, "%s.preopt", tree->io->out_lk_file);
		pfile = fopen(pfilestr3, "w");
		if(!pfile)
		{
			fprintf(stderr, "ERROR: cannot open output file: %s\n",
						pfilestr3);
			exit(1);
		}
		if(tree->io->print_site_lnl)
		{
			Print_Site_Lk(tree,pfile);
		}
		fclose(pfile);

		//
		// Print the pre-optimized sitecat lks
		//
		char sitecatlkfilestr[1000];
		strcpy(sitecatlkfilestr,tree->io->in_align_file);
		strcat(sitecatlkfilestr,"_phyml_sitecatlks");
	    if(tree->io->appebr_run_ID) { strcat(sitecatlkfilestr,"_"); strcat(sitecatlkfilestr,tree->io->run_id_string); }
	    strcat(sitecatlkfilestr,".txt.preopt");
	    FILE* sitecatfp = Openfile(sitecatlkfilestr,1);
		Print_Sitecat_Lks(tree, sitecatfp);
		fclose(sitecatfp);

		// Finally, return optimization algorithm to STA.
		tree->io->opt_algorithm = STA;
	}



	phydbl temp = tree->io->sa_temp_start;
	phydbl tempmult = exp(log(tree->io->sa_temp_end/tree->io->sa_temp_start)/(((phydbl)tree->io->sa_num_anneal_stages) - 1.0));
	if(tree->io->sa_scale_temp == 1)
	{
		anneal.t_scalar = Scale_T_Factor(tree);
		if (anneal.t_scalar < 0.1) // VHS: a hack?
		{
			anneal.t_scalar = 0.1;
		}
	}
	else{
		anneal.t_scalar = 1.0; //if positive, user input a value, don't estimate.
		Lk(tree);
	}


	//
	t_tree *best_tree = best_tree = Make_Tree_From_Scratch(tree->n_otu,tree->data);
	best_tree->mod  = Copy_Model(tree->mod);
#ifndef ENABLE_MBL
	phydbl *best_bl = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
#else
	phydbl *best_bl = (phydbl *)mCalloc((2*tree->n_otu-3)*n_l,sizeof(phydbl));
#endif
	Record_Br_Len(best_bl,tree);
	Copy_Tree(tree, best_tree);

	//
	t_tree *last_tree = Make_Tree_From_Scratch(tree->n_otu,tree->data);
	last_tree->mod = Copy_Model(tree->mod);
#ifndef ENABLE_MBL
	phydbl *last_bl = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
#else
	phydbl *last_bl = (phydbl *)mCalloc((2*tree->n_otu-3)*n_l,sizeof(phydbl));
#endif
	Record_Br_Len(last_bl,tree);
	Copy_Tree(tree, last_tree);

	time_t start = time(NULL);
	time_t now;

	Lk(tree);

	phydbl lnL_best = tree->c_lnL;
	phydbl lnL_current = tree->c_lnL;
	phydbl lnL_proposed = tree->c_lnL;
	phydbl acc_prob;
	phydbl r;
	int itemp,iter;
	int total_steps_tried = 0;
	int total_steps_accepted = 0;

	phydbl temp_of_best = 0; // VHS: for optimization purposes, let's record the temperature at which we discovered the best tree.
	phydbl iter_of_best = 0;


	for(itemp = 0; itemp < tree->io->sa_num_anneal_stages; itemp++)
	{
		int steps_tried = 0.0;
		int steps_accepted = 0.0;

		// Return to our best tree at the start of each temperature stage.
		Copy_Tree(best_tree,tree);
		Record_Model(best_tree->mod, tree->mod);
		Restore_Br_Len(best_bl, tree);

		Copy_Tree(tree,last_tree);
		Record_Model(best_tree->mod,last_tree->mod);
		Restore_Br_Len(best_bl, last_tree);

		lnL_current = lnL_best;

		for(iter = 0; iter < tree->io->sa_iters_per_stage; iter++)
		{
			steps_tried++;
			//get our next proposition.
			Get_Neighbor_Proposition(tree,temp);
			lnL_proposed = tree->c_lnL;

			now = time(NULL);
			phydbl save_lnL_current = (phydbl)lnL_current;

			if(lnL_proposed > lnL_best)
			{
				//save this tree into best_tree
				Copy_Tree(tree,best_tree); // write tree into best_tree
				Record_Model(tree->mod, best_tree->mod); // write tree->mod into best_tree->mod
				Record_Br_Len(best_bl, tree);

				lnL_best = lnL_proposed;
				temp_of_best = temp;
				iter_of_best = iter;

				// This temperature is yielding good results: let's stay here
				// longer, thus increasing our chances of reaching the best ground
				// state at this temperature.
				iter -= tree->io->sa_set_back;
				if(iter < 0) iter = 0; //make sure not set back into negative...
				acc_prob = 1.0;
				r = 0.0;
			}
			else if(lnL_proposed > lnL_current){
				iter -= tree->io->sa_set_back;
				if(iter < 0) iter = 0; //make sure not set back into negative...
				acc_prob = 1.0;
				r = 0.0;
			}else{
				acc_prob = Boltzmann_P_TA(lnL_current, lnL_proposed, temp);
				r = gsl_rng_uniform(anneal.rng);
			}

			if(acc_prob >= r){
				//save the current tree
				if(anneal.no_change == 0){
					Copy_Tree(tree,last_tree);
					Record_Model(tree->mod, last_tree->mod);
					Record_Br_Len(last_bl, last_tree);
				}
				steps_accepted++;
				lnL_current = tree->c_lnL;

				//PhyML_Printf("JSJ: Proposed Likelihood Accepted, acc_prob = %lf\n",acc_prob);

			}else{
				//restore the current tree to be our last tree
				if(anneal.no_change == 0){
					Copy_Tree(last_tree,tree);
					Record_Model(last_tree->mod, tree->mod);
					Restore_Br_Len(last_bl, tree);
				}
			}

			// calculate acceptance ratio
			phydbl acc_ratio = (phydbl)steps_accepted / (phydbl)steps_tried;

			PhyML_Printf("[%ld sec] [temp. = %f] [current lnL = %.4f] [proposed lnL = %.4f] [p = %.3f] [accepted/proposed = %.3f]\n", (long)difftime(now,start), temp, save_lnL_current, lnL_proposed, acc_prob, acc_ratio);

		}//end inner for loop.
		temp *= tempmult;

		total_steps_tried += steps_tried;
		total_steps_accepted += steps_accepted;


		// check for early convergence
		// i.e., there was no lnL change during the last temperature stage
		// the second half of this boolean ensures that we only stop if we've accomplished at least 1/3 of the cooling schedule
		if ( (tree->io->sa_stop_early == 1) && (itemp > ((phydbl)tree->io->sa_num_anneal_stages)/3.0) )
		{
		  if (steps_accepted == 0)
		    {
		      PhyML_Printf("\n. STA is stopping early because I think we reached convergence.");
		      break;
		    }
		}
	}
	if(verbose){
		PhyML_Printf("Annealing finished.\n");
		PhyML_Printf("In temperature range [%f, %f], the best tree was found at %f, iteration %f\n", temp, tree->io->sa_temp_start, temp_of_best, iter_of_best);
	}

	// Our final answer will the best_tree
	Copy_Tree(best_tree,tree);
	Record_Model(best_tree->mod, tree->mod);
	Restore_Br_Len(best_bl, tree);

	Free_Model(best_tree->mod);
	Free_Tree(best_tree);
	Free_Model(last_tree->mod);
	Free_Tree(last_tree);
	Free(best_bl);
	Free(last_bl);
	Free_Anneal(); //local function that frees the gsl_rng;
	Lk(tree);

#ifdef TIME
	clock_t time_STAstop = clock();
	time_STA += time_STAstop - time_STAstart;
#endif

	return tree->c_lnL;
}

#endif
